/*
 *  MWBoundaries.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 10-06-10.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include <cmath>
#include "MW.h"

void MW::Boundaries(Scalar *iPhi, Scalar *iPsi, Vector *iE, Vector *iB) {
	double RT2 = sqrt(2);
	double RT3 = sqrt(3);

	Scalar Phi = *iPhi;
	Scalar Psi = *iPsi;
	Vector E = *iE;
	Vector B = *iB;

	Scalar Phi2 = Phi;
	Scalar Psi2 = Psi;
	Vector E2 = E;
	Vector B2 = B;

	if ( not O.p(0) ) FOR(j,O[1]) FOR(k,O[2]) { // x=0,1 plane
		Phi2(0,j,k) = 0.5*(Phi(0,j,k)+E(0)(0,j,k));
		Psi2(0,j,k) = 0.5*(Psi(0,j,k)+B(0)(0,j,k));
		E2(0)(0,j,k) = 0.5*(E(0)(0,j,k)+Phi(0,j,k));
		E2(1)(0,j,k) = 0.5*(E(1)(0,j,k)-B(2)(0,j,k));
		E2(2)(0,j,k) = 0.5*(E(2)(0,j,k)+B(1)(0,j,k));
		B2(0)(0,j,k) = 0.5*(B(0)(0,j,k)+Psi(0,j,k));
		B2(1)(0,j,k) = 0.5*(B(1)(0,j,k)+E(2)(0,j,k));
		B2(2)(0,j,k) = 0.5*(B(2)(0,j,k)-E(1)(0,j,k));
		Phi2(-1,j,k) = 0.5*(Phi(-1,j,k)-E(0)(-1,j,k));
		Psi2(-1,j,k) = 0.5*(Psi(-1,j,k)-B(0)(-1,j,k));
		E2(0)(-1,j,k) = 0.5*(E(0)(-1,j,k)-Phi(-1,j,k));
		E2(1)(-1,j,k) = 0.5*(E(1)(-1,j,k)+B(2)(-1,j,k));
		E2(2)(-1,j,k) = 0.5*(E(2)(-1,j,k)-B(1)(-1,j,k));
		B2(0)(-1,j,k) = 0.5*(B(0)(-1,j,k)-Psi(-1,j,k));
		B2(1)(-1,j,k) = 0.5*(B(1)(-1,j,k)-E(2)(-1,j,k));
		B2(2)(-1,j,k) = 0.5*(B(2)(-1,j,k)+E(1)(-1,j,k));
	}
	if ( not O.p(1) ) FOR(i,O[0]) FOR(k,O[2]) { // y=0,1 plane
		Phi2(i,0,k) = 0.5*(Phi(i,0,k)+E(1)(i,0,k));
		Psi2(i,0,k) = 0.5*(Psi(i,0,k)+B(1)(i,0,k));
		E2(0)(i,0,k) = 0.5*(E(0)(i,0,k)+B(2)(i,0,k));
		E2(1)(i,0,k) = 0.5*(E(1)(i,0,k)+Phi(i,0,k));
		E2(2)(i,0,k) = 0.5*(E(2)(i,0,k)-B(0)(i,0,k));
		B2(0)(i,0,k) = 0.5*(B(0)(i,0,k)-E(2)(i,0,k));
		B2(1)(i,0,k) = 0.5*(B(1)(i,0,k)+Psi(i,0,k));
		B2(2)(i,0,k) = 0.5*(B(2)(i,0,k)+E(0)(i,0,k));
		Phi2(i,-1,k) = 0.5*(Phi(i,-1,k)-E(1)(i,-1,k));
		Psi2(i,-1,k) = 0.5*(Psi(i,-1,k)-B(1)(i,-1,k));
		E2(0)(i,-1,k) = 0.5*(E(0)(i,-1,k)-B(2)(i,-1,k));
		E2(1)(i,-1,k) = 0.5*(E(1)(i,-1,k)-Phi(i,-1,k));
		E2(2)(i,-1,k) = 0.5*(E(2)(i,-1,k)+B(0)(i,-1,k));
		B2(0)(i,-1,k) = 0.5*(B(0)(i,-1,k)+E(2)(i,-1,k));
		B2(1)(i,-1,k) = 0.5*(B(1)(i,-1,k)-Psi(i,-1,k));
		B2(2)(i,-1,k) = 0.5*(B(2)(i,-1,k)-E(0)(i,-1,k));
	}
	if ( not O.p(2) ) FOR(i,O[0]) FOR(j,O[1]) { // z=0,1 plane
		Phi2(i,j,0) = 0.5*(Phi(i,j,0)+E(2)(i,j,0));
		Psi2(i,j,0) = 0.5*(Psi(i,j,0)+B(2)(i,j,0));
		E2(0)(i,j,0) = 0.5*(E(0)(i,j,0)-B(1)(i,j,0));
		E2(1)(i,j,0) = 0.5*(E(1)(i,j,0)+B(0)(i,j,0));
		E2(2)(i,j,0) = 0.5*(E(2)(i,j,0)+Phi(i,j,0));
		B2(0)(i,j,0) = 0.5*(B(0)(i,j,0)+E(1)(i,j,0));
		B2(1)(i,j,0) = 0.5*(B(1)(i,j,0)-E(0)(i,j,0));
		B2(2)(i,j,0) = 0.5*(B(2)(i,j,0)+Psi(i,j,0));
		Phi2(i,j,-1) = 0.5*(Phi(i,j,-1)-E(2)(i,j,-1));
		Psi2(i,j,-1) = 0.5*(Psi(i,j,-1)-B(2)(i,j,-1));
		E2(0)(i,j,-1) = 0.5*(E(0)(i,j,-1)+B(1)(i,j,-1));
		E2(1)(i,j,-1) = 0.5*(E(1)(i,j,-1)-B(0)(i,j,-1));
		E2(2)(i,j,-1) = 0.5*(E(2)(i,j,-1)-Phi(i,j,-1));
		B2(0)(i,j,-1) = 0.5*(B(0)(i,j,-1)-E(1)(i,j,-1));
		B2(1)(i,j,-1) = 0.5*(B(1)(i,j,-1)+E(0)(i,j,-1));
		B2(2)(i,j,-1) = 0.5*(B(2)(i,j,-1)-Psi(i,j,-1));
	}

	if ( not O.p(0) and not O.p(1) ) FOR(k,O[2]) { // x,y = constant lines
		Phi2(0,0,k) = 0.5*(Phi(0,0,k)+(E(0)(0,0,k)+E(1)(0,0,k))/RT2);
		Psi2(0,0,k) = 0.5*(Psi(0,0,k)+(B(0)(0,0,k)+B(1)(0,0,k))/RT2);
		E2(0)(0,0,k) = 0.5*(E(0)(0,0,k)+(Phi(0,0,k)+B(2)(0,0,k))/RT2);
		E2(1)(0,0,k) = 0.5*(E(1)(0,0,k)+(Phi(0,0,k)-B(2)(0,0,k))/RT2);
		E2(2)(0,0,k) = 0.5*(E(2)(0,0,k)+(B(1)(0,0,k)-B(0)(0,0,k))/RT2);
		B2(0)(0,0,k) = 0.5*(B(0)(0,0,k)+(Psi(0,0,k)-E(2)(0,0,k))/RT2);
		B2(1)(0,0,k) = 0.5*(B(1)(0,0,k)+(Psi(0,0,k)+E(2)(0,0,k))/RT2);
		B2(2)(0,0,k) = 0.5*(B(2)(0,0,k)+(E(0)(0,0,k)-E(1)(0,0,k))/RT2);
		Phi2(-1,-1,k) = 0.5*(Phi(-1,-1,k)-(E(0)(-1,-1,k)+E(1)(-1,-1,k))/RT2);
		Psi2(-1,-1,k) = 0.5*(Psi(-1,-1,k)-(B(0)(-1,-1,k)+B(1)(-1,-1,k))/RT2);
		E2(0)(-1,-1,k) = 0.5*(E(0)(-1,-1,k)-(Phi(-1,-1,k)+B(2)(-1,-1,k))/RT2);
		E2(1)(-1,-1,k) = 0.5*(E(1)(-1,-1,k)-(Phi(-1,-1,k)-B(2)(-1,-1,k))/RT2);
		E2(2)(-1,-1,k) = 0.5*(E(2)(-1,-1,k)+(B(0)(-1,-1,k)-B(1)(-1,-1,k))/RT2);
		B2(0)(-1,-1,k) = 0.5*(B(0)(-1,-1,k)-(Psi(-1,-1,k)-E(2)(-1,-1,k))/RT2);
		B2(1)(-1,-1,k) = 0.5*(B(1)(-1,-1,k)-(Psi(-1,-1,k)+E(2)(-1,-1,k))/RT2);
		B2(2)(-1,-1,k) = 0.5*(B(2)(-1,-1,k)+(E(1)(-1,-1,k)-E(0)(-1,-1,k))/RT2);
		Phi2(-1,0,k) = 0.5*(Phi(-1,0,k)+(E(1)(-1,0,k)-E(0)(-1,0,k))/RT2);
		Psi2(-1,0,k) = 0.5*(Psi(-1,0,k)+(B(1)(-1,0,k)-B(0)(-1,0,k))/RT2);
		E2(0)(-1,0,k) = 0.5*(E(0)(-1,0,k)-(Phi(-1,0,k)-B(2)(-1,0,k))/RT2);
		E2(1)(-1,0,k) = 0.5*(E(1)(-1,0,k)+(Phi(-1,0,k)+B(2)(-1,0,k))/RT2);
		E2(2)(-1,0,k) = 0.5*(E(2)(-1,0,k)-(B(0)(-1,0,k)+B(1)(-1,0,k))/RT2);
		B2(0)(-1,0,k) = 0.5*(B(0)(-1,0,k)-(Psi(-1,0,k)+E(2)(-1,0,k))/RT2);
		B2(1)(-1,0,k) = 0.5*(B(1)(-1,0,k)+(Psi(-1,0,k)-E(2)(-1,0,k))/RT2);
		B2(2)(-1,0,k) = 0.5*(B(2)(-1,0,k)+(E(0)(-1,0,k)+E(1)(-1,0,k))/RT2);
		Phi2(0,-1,k) = 0.5*(Phi(0,-1,k)+(E(0)(0,-1,k)-E(1)(0,-1,k))/RT2);
		Psi2(0,-1,k) = 0.5*(Psi(0,-1,k)+(B(0)(0,-1,k)-B(1)(0,-1,k))/RT2);
		E2(0)(0,-1,k) = 0.5*(E(0)(0,-1,k)+(Phi(0,-1,k)-B(2)(0,-1,k))/RT2);
		E2(1)(0,-1,k) = 0.5*(E(1)(0,-1,k)-(Phi(0,-1,k)+B(2)(0,-1,k))/RT2);
		E2(2)(0,-1,k) = 0.5*(E(2)(0,-1,k)+(B(0)(0,-1,k)+B(1)(0,-1,k))/RT2);
		B2(0)(0,-1,k) = 0.5*(B(0)(0,-1,k)+(Psi(0,-1,k)+E(2)(0,-1,k))/RT2);
		B2(1)(0,-1,k) = 0.5*(B(1)(0,-1,k)-(Psi(0,-1,k)-E(2)(0,-1,k))/RT2);
		B2(2)(0,-1,k) = 0.5*(B(2)(0,-1,k)-(E(0)(0,-1,k)+E(1)(0,-1,k))/RT2);
	}

	if ( not O.p(1) and not O.p(2) ) FOR(i,O[0]) { // y,z = constant lines
		Phi2(i,0,0) = 0.5*(Phi(i,0,0)+(E(1)(i,0,0)+E(2)(i,0,0))/RT2);
		Psi2(i,0,0) = 0.5*(Psi(i,0,0)+(B(1)(i,0,0)+B(2)(i,0,0))/RT2);
		E2(0)(i,0,0) = 0.5*(E(0)(i,0,0)+(B(2)(i,0,0)-B(1)(i,0,0))/RT2);
		E2(1)(i,0,0) = 0.5*(E(1)(i,0,0)+(Phi(i,0,0)+B(0)(i,0,0))/RT2);
		E2(2)(i,0,0) = 0.5*(E(2)(i,0,0)+(Phi(i,0,0)-B(0)(i,0,0))/RT2);
		B2(0)(i,0,0) = 0.5*(B(0)(i,0,0)+(E(1)(i,0,0)-E(2)(i,0,0))/RT2);
		B2(1)(i,0,0) = 0.5*(B(1)(i,0,0)+(Psi(i,0,0)-E(0)(i,0,0))/RT2);
		B2(2)(i,0,0) = 0.5*(B(2)(i,0,0)+(Psi(i,0,0)+E(0)(i,0,0))/RT2);
		Phi2(i,-1,-1) = 0.5*(Phi(i,-1,-1)-(E(1)(i,-1,-1)+E(2)(i,-1,-1))/RT2);
		Psi2(i,-1,-1) = 0.5*(Psi(i,-1,-1)-(B(1)(i,-1,-1)+B(2)(i,-1,-1))/RT2);
		E2(0)(i,-1,-1) = 0.5*(E(0)(i,-1,-1)+(B(1)(i,-1,-1)-B(2)(i,-1,-1))/RT2);
		E2(1)(i,-1,-1) = 0.5*(E(1)(i,-1,-1)-(Phi(i,-1,-1)+B(0)(i,-1,-1))/RT2);
		E2(2)(i,-1,-1) = 0.5*(E(2)(i,-1,-1)-(Phi(i,-1,-1)-B(0)(i,-1,-1))/RT2);
		B2(0)(i,-1,-1) = 0.5*(B(0)(i,-1,-1)+(E(2)(i,-1,-1)-E(1)(i,-1,-1))/RT2);
		B2(1)(i,-1,-1) = 0.5*(B(1)(i,-1,-1)-(Psi(i,-1,-1)-E(0)(i,-1,-1))/RT2);
		B2(2)(i,-1,-1) = 0.5*(B(2)(i,-1,-1)-(Psi(i,-1,-1)+E(0)(i,-1,-1))/RT2);
		Phi2(i,-1,0) = 0.5*(Phi(i,-1,0)+(E(2)(i,-1,0)-E(1)(i,-1,0))/RT2);
		Psi2(i,-1,0) = 0.5*(Psi(i,-1,0)+(B(2)(i,-1,0)-B(1)(i,-1,0))/RT2);
		E2(0)(i,-1,0) = 0.5*(E(0)(i,-1,0)-(B(1)(i,-1,0)+B(2)(i,-1,0))/RT2);
		E2(1)(i,-1,0) = 0.5*(E(1)(i,-1,0)-(Phi(i,-1,0)-B(0)(i,-1,0))/RT2);
		E2(2)(i,-1,0) = 0.5*(E(2)(i,-1,0)+(Phi(i,-1,0)+B(0)(i,-1,0))/RT2);
		B2(0)(i,-1,0) = 0.5*(B(0)(i,-1,0)+(E(1)(i,-1,0)+E(2)(i,-1,0))/RT2);
		B2(1)(i,-1,0) = 0.5*(B(1)(i,-1,0)-(Psi(i,-1,0)+E(0)(i,-1,0))/RT2);
		B2(2)(i,-1,0) = 0.5*(B(2)(i,-1,0)+(Psi(i,-1,0)-E(0)(i,-1,0))/RT2);
		Phi2(i,0,-1) = 0.5*(Phi(i,0,-1)+(E(1)(i,0,-1)-E(2)(i,0,-1))/RT2);
		Psi2(i,0,-1) = 0.5*(Psi(i,0,-1)+(B(1)(i,0,-1)-B(2)(i,0,-1))/RT2);
		E2(0)(i,0,-1) = 0.5*(E(0)(i,0,-1)+(B(1)(i,0,-1)+B(2)(i,0,-1))/RT2);
		E2(1)(i,0,-1) = 0.5*(E(1)(i,0,-1)+(Phi(i,0,-1)-B(0)(i,0,-1))/RT2);
		E2(2)(i,0,-1) = 0.5*(E(2)(i,0,-1)-(Phi(i,0,-1)+B(0)(i,0,-1))/RT2);
		B2(0)(i,0,-1) = 0.5*(B(0)(i,0,-1)-(E(1)(i,0,-1)+E(2)(i,0,-1))/RT2);
		B2(1)(i,0,-1) = 0.5*(B(1)(i,0,-1)+(Psi(i,0,-1)+E(0)(i,0,-1))/RT2);
		B2(2)(i,0,-1) = 0.5*(B(2)(i,0,-1)-(Psi(i,0,-1)-E(0)(i,0,-1))/RT2);
	}

	if ( not O.p(0) and not O.p(2) ) FOR(j,O[1]) { // x,z = constant lines
		Phi2(0,j,0) = 0.5*(Phi(0,j,0)+(E(2)(0,j,0)+E(0)(0,j,0))/RT2);
		Psi2(0,j,0) = 0.5*(Psi(0,j,0)+(B(2)(0,j,0)+B(0)(0,j,0))/RT2);
		E2(0)(0,j,0) = 0.5*(E(0)(0,j,0)+(Phi(0,j,0)-B(1)(0,j,0))/RT2);
		E2(1)(0,j,0) = 0.5*(E(1)(0,j,0)+(B(0)(0,j,0)-B(2)(0,j,0))/RT2);
		E2(2)(0,j,0) = 0.5*(E(2)(0,j,0)+(Phi(0,j,0)+B(1)(0,j,0))/RT2);
		B2(0)(0,j,0) = 0.5*(B(0)(0,j,0)+(Psi(0,j,0)+E(1)(0,j,0))/RT2);
		B2(1)(0,j,0) = 0.5*(B(1)(0,j,0)+(E(2)(0,j,0)-E(0)(0,j,0))/RT2);
		B2(2)(0,j,0) = 0.5*(B(2)(0,j,0)+(Psi(0,j,0)-E(1)(0,j,0))/RT2);
		Phi2(-1,j,-1) = 0.5*(Phi(-1,j,-1)-(E(2)(-1,j,-1)+E(0)(-1,j,-1))/RT2);
		Psi2(-1,j,-1) = 0.5*(Psi(-1,j,-1)-(B(2)(-1,j,-1)+B(0)(-1,j,-1))/RT2);
		E2(0)(-1,j,-1) = 0.5*(E(0)(-1,j,-1)-(Phi(-1,j,-1)-B(1)(-1,j,-1))/RT2);
		E2(1)(-1,j,-1) = 0.5*(E(1)(-1,j,-1)+(B(2)(-1,j,-1)-B(0)(-1,j,-1))/RT2);
		E2(2)(-1,j,-1) = 0.5*(E(2)(-1,j,-1)-(Phi(-1,j,-1)+B(1)(-1,j,-1))/RT2);
		B2(0)(-1,j,-1) = 0.5*(B(0)(-1,j,-1)-(Psi(-1,j,-1)+E(1)(-1,j,-1))/RT2);
		B2(1)(-1,j,-1) = 0.5*(B(1)(-1,j,-1)+(E(0)(-1,j,-1)-E(2)(-1,j,-1))/RT2);
		B2(2)(-1,j,-1) = 0.5*(B(2)(-1,j,-1)-(Psi(-1,j,-1)-E(1)(-1,j,-1))/RT2);
		Phi2(0,j,-1) = 0.5*(Phi(0,j,-1)+(E(0)(0,j,-1)-E(2)(0,j,-1))/RT2);
		Psi2(0,j,-1) = 0.5*(Psi(0,j,-1)+(B(0)(0,j,-1)-B(2)(0,j,-1))/RT2);
		E2(0)(0,j,-1) = 0.5*(E(0)(0,j,-1)+(Phi(0,j,-1)+B(1)(0,j,-1))/RT2);
		E2(1)(0,j,-1) = 0.5*(E(1)(0,j,-1)-(B(2)(0,j,-1)+B(0)(0,j,-1))/RT2);
		E2(2)(0,j,-1) = 0.5*(E(2)(0,j,-1)-(Phi(0,j,-1)-B(1)(0,j,-1))/RT2);
		B2(0)(0,j,-1) = 0.5*(B(0)(0,j,-1)+(Psi(0,j,-1)-E(1)(0,j,-1))/RT2);
		B2(1)(0,j,-1) = 0.5*(B(1)(0,j,-1)+(E(2)(0,j,-1)+E(0)(0,j,-1))/RT2);
		B2(2)(0,j,-1) = 0.5*(B(2)(0,j,-1)-(Psi(0,j,-1)+E(1)(0,j,-1))/RT2);
		Phi2(-1,j,0) = 0.5*(Phi(-1,j,0)+(E(2)(-1,j,0)-E(0)(-1,j,0))/RT2);
		Psi2(-1,j,0) = 0.5*(Psi(-1,j,0)+(B(2)(-1,j,0)-B(0)(-1,j,0))/RT2);
		E2(0)(-1,j,0) = 0.5*(E(0)(-1,j,0)-(Phi(-1,j,0)+B(1)(-1,j,0))/RT2);
		E2(1)(-1,j,0) = 0.5*(E(1)(-1,j,0)+(B(2)(-1,j,0)+B(0)(-1,j,0))/RT2);
		E2(2)(-1,j,0) = 0.5*(E(2)(-1,j,0)+(Phi(-1,j,0)-B(1)(-1,j,0))/RT2);
		B2(0)(-1,j,0) = 0.5*(B(0)(-1,j,0)-(Psi(-1,j,0)-E(1)(-1,j,0))/RT2);
		B2(1)(-1,j,0) = 0.5*(B(1)(-1,j,0)-(E(2)(-1,j,0)+E(0)(-1,j,0))/RT2);
		B2(2)(-1,j,0) = 0.5*(B(2)(-1,j,0)+(Psi(-1,j,0)+E(1)(-1,j,0))/RT2);
	}

	if ( not O.p(0) and not O.p(1) and not O.p(2) ) {
		{ // point 0,0,0
			Phi2(0,0,0) = 0.5*(Phi(0,0,0)+(E(0)(0,0,0)+E(1)(0,0,0)+E(2)(0,0,0))/RT3);
			Psi2(0,0,0) = 0.5*(Psi(0,0,0)+(B(0)(0,0,0)+B(1)(0,0,0)+B(2)(0,0,0))/RT3);
			E2(0)(0,0,0) = 0.5*(E(0)(0,0,0)+(Phi(0,0,0)-B(1)(0,0,0)+B(2)(0,0,0))/RT3);
			E2(1)(0,0,0) = 0.5*(E(1)(0,0,0)+(Phi(0,0,0)+B(0)(0,0,0)-B(2)(0,0,0))/RT3);
			E2(2)(0,0,0) = 0.5*(E(2)(0,0,0)+(Phi(0,0,0)-B(0)(0,0,0)+B(1)(0,0,0))/RT3);
			B2(0)(0,0,0) = 0.5*(B(0)(0,0,0)+(Psi(0,0,0)+E(1)(0,0,0)-E(2)(0,0,0))/RT3);
			B2(1)(0,0,0) = 0.5*(B(1)(0,0,0)+(Psi(0,0,0)-E(0)(0,0,0)+E(2)(0,0,0))/RT3);
			B2(2)(0,0,0) = 0.5*(B(2)(0,0,0)+(Psi(0,0,0)+E(0)(0,0,0)-E(1)(0,0,0))/RT3);
		}
		{ // point 0,0,1
			Phi2(0,0,-1) = 0.5*(Phi(0,0,-1)+(E(0)(0,0,-1)+E(1)(0,0,-1)-E(2)(0,0,-1))/RT3);
			Psi2(0,0,-1) = 0.5*(Psi(0,0,-1)+(B(0)(0,0,-1)+B(1)(0,0,-1)-B(2)(0,0,-1))/RT3);
			E2(0)(0,0,-1) = 0.5*(E(0)(0,0,-1)+(Phi(0,0,-1)+B(1)(0,0,-1)+B(2)(0,0,-1))/RT3);
			E2(1)(0,0,-1) = 0.5*(E(1)(0,0,-1)+(Phi(0,0,-1)-B(0)(0,0,-1)-B(2)(0,0,-1))/RT3);
			E2(2)(0,0,-1) = 0.5*(E(2)(0,0,-1)-(Phi(0,0,-1)+B(0)(0,0,-1)-B(1)(0,0,-1))/RT3);
			B2(0)(0,0,-1) = 0.5*(B(0)(0,0,-1)+(Psi(0,0,-1)-E(1)(0,0,-1)-E(2)(0,0,-1))/RT3);
			B2(1)(0,0,-1) = 0.5*(B(1)(0,0,-1)+(Psi(0,0,-1)+E(0)(0,0,-1)+E(2)(0,0,-1))/RT3);
			B2(2)(0,0,-1) = 0.5*(B(2)(0,0,-1)-(Psi(0,0,-1)-E(0)(0,0,-1)+E(1)(0,0,-1))/RT3);
		}
		{ // point 0,1,0
			Phi2(0,-1,0) = 0.5*(Phi(0,-1,0)+(E(0)(0,-1,0)-E(1)(0,-1,0)+E(2)(0,-1,0))/RT3);
			Psi2(0,-1,0) = 0.5*(Psi(0,-1,0)+(B(0)(0,-1,0)-B(1)(0,-1,0)+B(2)(0,-1,0))/RT3);
			E2(0)(0,-1,0) = 0.5*(E(0)(0,-1,0)+(Phi(0,-1,0)-B(1)(0,-1,0)-B(2)(0,-1,0))/RT3);
			E2(1)(0,-1,0) = 0.5*(E(1)(0,-1,0)-(Phi(0,-1,0)-B(0)(0,-1,0)+B(2)(0,-1,0))/RT3);
			E2(2)(0,-1,0) = 0.5*(E(2)(0,-1,0)+(Phi(0,-1,0)+B(0)(0,-1,0)+B(1)(0,-1,0))/RT3);
			B2(0)(0,-1,0) = 0.5*(B(0)(0,-1,0)+(Psi(0,-1,0)+E(1)(0,-1,0)+E(2)(0,-1,0))/RT3);
			B2(1)(0,-1,0) = 0.5*(B(1)(0,-1,0)-(Psi(0,-1,0)+E(0)(0,-1,0)-E(2)(0,-1,0))/RT3);
			B2(2)(0,-1,0) = 0.5*(B(2)(0,-1,0)+(Psi(0,-1,0)-E(0)(0,-1,0)-E(1)(0,-1,0))/RT3);
		}
		{ // point 1,0,0
			Phi2(-1,0,0) = 0.5*(Phi(-1,0,0)+(-E(0)(-1,0,0)+E(1)(-1,0,0)+E(2)(-1,0,0))/RT3);
			Psi2(-1,0,0) = 0.5*(Psi(-1,0,0)+(-B(0)(-1,0,0)+B(1)(-1,0,0)+B(2)(-1,0,0))/RT3);
			E2(0)(-1,0,0) = 0.5*(E(0)(-1,0,0)-(Phi(-1,0,0)+B(1)(-1,0,0)-B(2)(-1,0,0))/RT3);
			E2(1)(-1,0,0) = 0.5*(E(1)(-1,0,0)+(Phi(-1,0,0)+B(0)(-1,0,0)+B(2)(-1,0,0))/RT3);
			E2(2)(-1,0,0) = 0.5*(E(2)(-1,0,0)+(Phi(-1,0,0)-B(0)(-1,0,0)-B(1)(-1,0,0))/RT3);
			B2(0)(-1,0,0) = 0.5*(B(0)(-1,0,0)-(Psi(-1,0,0)-E(1)(-1,0,0)+E(2)(-1,0,0))/RT3);
			B2(1)(-1,0,0) = 0.5*(B(1)(-1,0,0)+(Psi(-1,0,0)-E(0)(-1,0,0)-E(2)(-1,0,0))/RT3);
			B2(2)(-1,0,0) = 0.5*(B(2)(-1,0,0)+(Psi(-1,0,0)+E(0)(-1,0,0)+E(1)(-1,0,0))/RT3);
		}
		{ // point 1,1,0
			Phi2(-1,-1,0) = 0.5*(Phi(-1,-1,0)+(-E(0)(-1,-1,0)-E(1)(-1,-1,0)+E(2)(-1,-1,0))/RT3);
			Psi2(-1,-1,0) = 0.5*(Psi(-1,-1,0)+(-B(0)(-1,-1,0)-B(1)(-1,-1,0)+B(2)(-1,-1,0))/RT3);
			E2(0)(-1,-1,0) = 0.5*(E(0)(-1,-1,0)-(Phi(-1,-1,0)+B(1)(-1,-1,0)+B(2)(-1,-1,0))/RT3);
			E2(1)(-1,-1,0) = 0.5*(E(1)(-1,-1,0)-(Phi(-1,-1,0)-B(0)(-1,-1,0)-B(2)(-1,-1,0))/RT3);
			E2(2)(-1,-1,0) = 0.5*(E(2)(-1,-1,0)+(Phi(-1,-1,0)+B(0)(-1,-1,0)-B(1)(-1,-1,0))/RT3);
			B2(0)(-1,-1,0) = 0.5*(B(0)(-1,-1,0)-(Psi(-1,-1,0)-E(1)(-1,-1,0)-E(2)(-1,-1,0))/RT3);
			B2(1)(-1,-1,0) = 0.5*(B(1)(-1,-1,0)-(Psi(-1,-1,0)+E(0)(-1,-1,0)+E(2)(-1,-1,0))/RT3);
			B2(2)(-1,-1,0) = 0.5*(B(2)(-1,-1,0)+(Psi(-1,-1,0)-E(0)(-1,-1,0)+E(1)(-1,-1,0))/RT3);
		}
		{ // point 1,0,1
			Phi2(-1,0,-1) = 0.5*(Phi(-1,0,-1)+(-E(0)(-1,0,-1)+E(1)(-1,0,-1)-E(2)(-1,0,-1))/RT3);
			Psi2(-1,0,-1) = 0.5*(Psi(-1,0,-1)+(-B(0)(-1,0,-1)+B(1)(-1,0,-1)-B(2)(-1,0,-1))/RT3);
			E2(0)(-1,0,-1) = 0.5*(E(0)(-1,0,-1)-(Phi(-1,0,-1)-B(1)(-1,0,-1)-B(2)(-1,0,-1))/RT3);
			E2(1)(-1,0,-1) = 0.5*(E(1)(-1,0,-1)+(Phi(-1,0,-1)-B(0)(-1,0,-1)+B(2)(-1,0,-1))/RT3);
			E2(2)(-1,0,-1) = 0.5*(E(2)(-1,0,-1)-(Phi(-1,0,-1)+B(0)(-1,0,-1)+B(1)(-1,0,-1))/RT3);
			B2(0)(-1,0,-1) = 0.5*(B(0)(-1,0,-1)-(Psi(-1,0,-1)+E(1)(-1,0,-1)+E(2)(-1,0,-1))/RT3);
			B2(1)(-1,0,-1) = 0.5*(B(1)(-1,0,-1)+(Psi(-1,0,-1)+E(0)(-1,0,-1)-E(2)(-1,0,-1))/RT3);
			B2(2)(-1,0,-1) = 0.5*(B(2)(-1,0,-1)-(Psi(-1,0,-1)-E(0)(-1,0,-1)-E(1)(-1,0,-1))/RT3);
		}
		{ // point 0,1,1
			Phi2(0,-1,-1) = 0.5*(Phi(0,-1,-1)+(E(0)(0,-1,-1)-E(1)(0,-1,-1)-E(2)(0,-1,-1))/RT3);
			Psi2(0,-1,-1) = 0.5*(Psi(0,-1,-1)+(B(0)(0,-1,-1)-B(1)(0,-1,-1)-B(2)(0,-1,-1))/RT3);
			E2(0)(0,-1,-1) = 0.5*(E(0)(0,-1,-1)+(Phi(0,-1,-1)+B(1)(0,-1,-1)-B(2)(0,-1,-1))/RT3);
			E2(1)(0,-1,-1) = 0.5*(E(1)(0,-1,-1)-(Phi(0,-1,-1)+B(0)(0,-1,-1)+B(2)(0,-1,-1))/RT3);
			E2(2)(0,-1,-1) = 0.5*(E(2)(0,-1,-1)-(Phi(0,-1,-1)-B(0)(0,-1,-1)-B(1)(0,-1,-1))/RT3);
			B2(0)(0,-1,-1) = 0.5*(B(0)(0,-1,-1)+(Psi(0,-1,-1)-E(1)(0,-1,-1)+E(2)(0,-1,-1))/RT3);
			B2(1)(0,-1,-1) = 0.5*(B(1)(0,-1,-1)-(Psi(0,-1,-1)-E(0)(0,-1,-1)-E(2)(0,-1,-1))/RT3);
			B2(2)(0,-1,-1) = 0.5*(B(2)(0,-1,-1)-(Psi(0,-1,-1)+E(0)(0,-1,-1)+E(1)(0,-1,-1))/RT3);
		}
		{ // point 1,1,1
			Phi2(-1,-1,-1) = 0.5*(Phi(-1,-1,-1)-(E(0)(-1,-1,-1)+E(1)(-1,-1,-1)+E(2)(-1,-1,-1))/RT3);
			Psi2(-1,-1,-1) = 0.5*(Psi(-1,-1,-1)-(B(0)(-1,-1,-1)+B(1)(-1,-1,-1)+B(2)(-1,-1,-1))/RT3);
			E2(0)(-1,-1,-1) = 0.5*(E(0)(-1,-1,-1)-(Phi(-1,-1,-1)-B(1)(-1,-1,-1)+B(2)(-1,-1,-1))/RT3);
			E2(1)(-1,-1,-1) = 0.5*(E(1)(-1,-1,-1)-(Phi(-1,-1,-1)+B(0)(-1,-1,-1)-B(2)(-1,-1,-1))/RT3);
			E2(2)(-1,-1,-1) = 0.5*(E(2)(-1,-1,-1)-(Phi(-1,-1,-1)-B(0)(-1,-1,-1)+B(1)(-1,-1,-1))/RT3);
			B2(0)(-1,-1,-1) = 0.5*(B(0)(-1,-1,-1)-(Psi(-1,-1,-1)+E(1)(-1,-1,-1)-E(2)(-1,-1,-1))/RT3);
			B2(1)(-1,-1,-1) = 0.5*(B(1)(-1,-1,-1)-(Psi(-1,-1,-1)-E(0)(-1,-1,-1)+E(2)(-1,-1,-1))/RT3);
			B2(2)(-1,-1,-1) = 0.5*(B(2)(-1,-1,-1)-(Psi(-1,-1,-1)+E(0)(-1,-1,-1)-E(1)(-1,-1,-1))/RT3);
		}
	}

	if ( O.p(0) ) FOR(j,O[1]) FOR(k,O[2]) {
		Phi2(0,j,k) = Phi2(-1,j,k) = 0.5*(Phi(0,j,k)+Phi(-1,j,k));
		Psi2(0,j,k) = Psi2(-1,j,k) = 0.5*(Psi(0,j,k)+Psi(-1,j,k));
		FORa E2(a)(0,j,k) = E2(a)(-1,j,k) = 0.5*(E(a)(0,j,k)+E(a)(-1,j,k));
		FORa B2(a)(0,j,k) = B2(a)(-1,j,k) = 0.5*(B(a)(0,j,k)+B(a)(-1,j,k));
	}
	if ( O.p(1) ) FOR(i,O[0]) FOR(k,O[2]) {
		Phi2(i,0,k) = Phi2(i,-1,k) = 0.5*(Phi(i,0,k)+Phi(i,-1,k));
		Psi2(i,0,k) = Psi2(i,-1,k) = 0.5*(Psi(i,0,k)+Psi(i,-1,k));
		FORa E2(a)(i,0,k) = E2(a)(i,-1,k) = 0.5*(E(a)(i,0,k)+E(a)(i,-1,k));
		FORa B2(a)(i,0,k) = B2(a)(i,-1,k) = 0.5*(B(a)(i,0,k)+B(a)(i,-1,k));
	}
	if ( O.p(2) ) FOR(i,O[0]) FOR(j,O[1]) {
		Phi2(i,j,0) = Phi2(i,j,-1) = 0.5*(Phi(i,j,0)+Phi(i,j,-1));
		Psi2(i,j,0) = Psi2(i,j,-1) = 0.5*(Psi(i,j,0)+Psi(i,j,-1));
		FORa E2(a)(i,j,0) = E2(a)(i,j,-1) = 0.5*(E(a)(i,j,0)+E(a)(i,j,-1));
		FORa B2(a)(i,j,0) = B2(a)(i,j,-1) = 0.5*(B(a)(i,j,0)+B(a)(i,j,-1));
	}
	if ( O.p(0) and O.p(1) ) FOR(k,O[2]) {
		Phi2(0,0,k) = Phi2(0,-1,k) = Phi2(-1,0,k) = Phi2(-1,-1,k) = 0.25*(Phi(0,0,k)+Phi(0,-1,k)+Phi(-1,0,k)+Phi(-1,-1,k));
		Psi2(0,0,k) = Psi2(0,-1,k) = Psi2(-1,0,k) = Psi2(-1,-1,k) = 0.25*(Psi(0,0,k)+Psi(0,-1,k)+Psi(-1,0,k)+Psi(-1,-1,k));
		FORa E2(a)(0,0,k) = E2(a)(0,-1,k) = E2(a)(-1,0,k) = E2(a)(-1,-1,k) = 0.25*(E(a)(0,0,k)+E(a)(0,-1,k)+E(a)(-1,0,k)+E(a)(-1,-1,k));
		FORa B2(a)(0,0,k) = B2(a)(0,-1,k) = B2(a)(-1,0,k) = B2(a)(-1,-1,k) = 0.25*(B(a)(0,0,k)+B(a)(0,-1,k)+B(a)(-1,0,k)+B(a)(-1,-1,k));
	}
	if ( O.p(1) and O.p(2) ) FOR(i,O[0]) {
		Phi2(i,0,0) = Phi2(i,0,-1) = Phi2(i,-1,0) = Phi2(i,-1,-1) = 0.25*(Phi(i,0,0)+Phi(i,0,-1)+Phi(i,-1,0)+Phi(i,-1,-1));
		Psi2(i,0,0) = Psi2(i,0,-1) = Psi2(i,-1,0) = Psi2(i,-1,-1) = 0.25*(Psi(i,0,0)+Psi(i,0,-1)+Psi(i,-1,0)+Psi(i,-1,-1));
		FORa E2(a)(i,0,0) = E2(a)(i,0,-1) = E2(a)(i,-1,0) = E2(a)(i,-1,-1) = 0.25*(E(a)(i,0,0)+E(a)(i,0,-1)+E(a)(i,-1,0)+E(a)(i,-1,-1));
		FORa B2(a)(i,0,0) = B2(a)(i,0,-1) = B2(a)(i,-1,0) = B2(a)(i,-1,-1) = 0.25*(B(a)(i,0,0)+B(a)(i,0,-1)+B(a)(i,-1,0)+B(a)(i,-1,-1));
	}
	if ( O.p(0) and O.p(2) ) FOR(j,O[1]) {
		Phi2(0,j,0) = Phi2(0,j,-1) = Phi2(-1,j,0) = Phi2(-1,j,-1) = 0.25*(Phi(0,j,0)+Phi(0,j,-1)+Phi(-1,j,0)+Phi(-1,j,-1));
		Psi2(0,j,0) = Psi2(0,j,-1) = Psi2(-1,j,0) = Psi2(-1,j,-1) = 0.25*(Psi(0,j,0)+Psi(0,j,-1)+Psi(-1,j,0)+Psi(-1,j,-1));
		FORa E2(a)(0,j,0) = E2(a)(0,j,-1) = E2(a)(-1,j,0) = E2(a)(-1,j,-1) = 0.25*(E(a)(0,j,0)+E(a)(0,j,-1)+E(a)(-1,j,0)+E(a)(-1,j,-1));
		FORa B2(a)(0,j,0) = B2(a)(0,j,-1) = B2(a)(-1,j,0) = B2(a)(-1,j,-1) = 0.25*(B(a)(0,j,0)+B(a)(0,j,-1)+B(a)(-1,j,0)+B(a)(-1,j,-1));
	}
	if ( O.p(0) and O.p(1) and O.p(2) ) {
		Phi2(0,0,0) = Phi2(0,0,-1) = Phi2(0,-1,0) = Phi2(0,-1,-1) = Phi2(-1,0,0) = Phi2(-1,0,-1) = Phi2(-1,-1,0) = Phi2(-1,-1,-1) = 0.125*(Phi(0,0,0)+Phi(0,0,-1)+Phi(0,-1,0)+Phi(0,-1,-1)+Phi(-1,0,0)+Phi(-1,0,-1)+Phi(-1,-1,0)+Phi(-1,-1,-1));
		Psi2(0,0,0) = Psi2(0,0,-1) = Psi2(0,-1,0) = Psi2(0,-1,-1) = Psi2(-1,0,0) = Psi2(-1,0,-1) = Psi2(-1,-1,0) = Psi2(-1,-1,-1) = 0.125*(Psi(0,0,0)+Psi(0,0,-1)+Psi(0,-1,0)+Psi(0,-1,-1)+Psi(-1,0,0)+Psi(-1,0,-1)+Psi(-1,-1,0)+Psi(-1,-1,-1));
		FORa E2(a)(0,0,0) = E2(a)(0,0,-1) = E2(a)(0,-1,0) = E2(a)(0,-1,-1) = E2(a)(-1,0,0) = E2(a)(-1,0,-1) = E2(a)(-1,-1,0) = E2(a)(-1,-1,-1) = 0.125*(E(a)(0,0,0)+E(a)(0,0,-1)+E(a)(0,-1,0)+E(a)(0,-1,-1)+E(a)(-1,0,0)+E(a)(-1,0,-1)+E(a)(-1,-1,0)+E(a)(-1,-1,-1));
		FORa B2(a)(0,0,0) = B2(a)(0,0,-1) = B2(a)(0,-1,0) = B2(a)(0,-1,-1) = B2(a)(-1,0,0) = B2(a)(-1,0,-1) = B2(a)(-1,-1,0) = B2(a)(-1,-1,-1) = 0.125*(B(a)(0,0,0)+B(a)(0,0,-1)+B(a)(0,-1,0)+B(a)(0,-1,-1)+B(a)(-1,0,0)+B(a)(-1,0,-1)+B(a)(-1,-1,0)+B(a)(-1,-1,-1));
	}


	*iPhi = Phi2;
	*iPsi = Psi2;
	*iE = E2;
	*iB = B2;
}

